/** Fichier implémentant un ellipsoide.
 *
 * Ce fichier contient l'implémentation de la classe Ellipsoide.
 *
 * @file Ellipsoide.cpp
 * @author Johan "Solo" Luisier
 * @date 23 Juin 2011
 */
#include "Ellipsoide.hpp"

#include "Algebre/Solveur.hpp"

#include "Utils/Erreur.hpp"

#include <cmath>

#include <stdlib.h>

using namespace std;

/**  Constructeur par défaut.
 *
 * Ce constructeur crée un Ellipsoide centré à l'origine et de rayon 1. selon
 * Oy et Oz, de rayon 2. selon Ox.
 */
Ellipsoide::Ellipsoide()
    : Forme( "Ellipsoïde" ), Centre( 0., 0., 0. ), AxeX( 1., 0., 0. ),
      AxeY( 0., 1., 0. ), AxeZ( 0., 0., 1. ), RayonX( 2. ), RayonY( 1. ),
      RayonZ( 1. ), VersSphere( .5, 1., 1. )
{
    PseudoCentre = Centre;
    PseudoRayon = max( max( RayonX, RayonY ), RayonZ );
    MoteurResolution = new Solveur();

    VersLocal = Transformation( MatriceSO3(), Centre );
    VersGlobal = VersLocal.inverse();
}

/** Constructeur standard.
 *
 * Ce constructeur initialise le centre de l'Ellipsoide, ses axes
 * principaux et ses rayons. Un rayon négatif est automatiquement
 * transformé en valeur positive.
 *
 * @param[in] pt centre de l'ellipsoïde.
 * @param[in] x axe Ox de l'ellipsoïde.
 * @param[in] y axe Oy de l'ellipsoïde.
 * @param[in] rx rayon de l'ellipsoïde selon son axe Ox.
 * @param[in] ry rayon de l'ellipsoïde selon son axe Oy.
 * @param[in] rz rayon de l'ellipsoïde selon son axe Oz.
 */
Ellipsoide::Ellipsoide( const Coord& pt, const VecteurU& x, const VecteurU& y,
                        const double& rx, const double& ry, const double& rz )
    : Forme( "Ellipsoïde" ), Centre( pt ),
      AxeX( x ), AxeY( y ), AxeZ( x ^ y ),
      RayonX( abs( rx ) ), RayonY( abs( ry ) ), RayonZ( abs( rz ) ),
      VersSphere( 1. / rx, 1. / ry, 1. / rz )
{
    PseudoCentre = Centre;
    PseudoRayon =  max( max( rx, ry ), rz );
    MoteurResolution = new Solveur();

    if ( AxeZ.norme() < 1.e-14 )
    {
        cerr << "Les deux axes donnés sont collinéaires!" << endl;
        exit( 1 );
    }

    MatriceSO3 * test = 0;

    try
    {
        test = new MatriceSO3( AxeX, AxeY, AxeZ );
    }
    catch ( Erreur& e )
    {
        cerr << e.message << endl;
        delete test;
        exit( 1 );
    }

    VersLocal = Transformation( *test, Centre );
    VersGlobal = VersLocal.inverse();
    delete test;
}

/** Constructeur de copie, permet de dupliquer un Ellipsoide.
 *
 * @param[in] s objet à copier.
 */
Ellipsoide::Ellipsoide(const Ellipsoide &s)
    : Forme( "Ellipsoïde", s.VersGlobal ), Centre( s.Centre ),
      AxeX( s.AxeX ), AxeY( s.AxeY ), AxeZ( s.AxeZ ),
      RayonX( s.RayonX ), RayonY( s.RayonY ), RayonZ( s.RayonZ )
    
{
    PseudoCentre = Centre;
    MoteurResolution = new Solveur();
}

/** Destructeur, virtuel pour assurer une séquence de destruction correcte.
 *
 * Il faut libérer la mémoire allouée pour le solveur.
 */
Ellipsoide::~Ellipsoide()
{}

/** Accesseur pour RayonX.
 *
 * Permet de récupérer la valeur du rayon sur l'axe Ox.
 *
 * @return une référence constante sur RayonX.
 */
const double& Ellipsoide::rayonX() const
{
    return RayonX;
}

/** Accesseur pour RayonY.
 *
 * Permet de récupérer la valeur du rayon sur l'axe Oy.
 *
 * @return une référence constante sur RayonY.
 */
const double& Ellipsoide::rayonY() const
{
    return RayonY;
}

/** Accesseur pour RayonZ.
 *
 * Permet de récupérer la valeur du rayon sur l'axe Oz.
 *
 * @return une référence constante sur RayonZ.
 */
const double& Ellipsoide::rayonZ() const
{
    return RayonZ;
}

/** Accesseur pour Centre.
 *
 * Permet de récupérer le centre de l'Ellipsoide.
 *
 * @return une référence constante sur Centre.
 */
const Coord& Ellipsoide::centre() const
{
    return Centre;
}

/** Calcule l'intersection entre le rayon et l'ellispoïde formé par l'instance.
 *
 * Calcule l'intersection d'un rayon avec l'ellispoïde. Pour se faire on
 * procède à deux simplifications~: premièrement on passe dans le système de
 * coordonées locales (dans lequel l'elipsoïde est centré sur l'origine) et on
 * effectue ensuite une homothétie pour se ramener au cas d'une sphère unité.
 *
 * @param[in] ray instance de rayon avec lequel il faut trouver l'intersection.
 *
 * @retval NullePart si l'intersection n'existe pas, @retval inter les
 * coordonnées du point d'intersection.
 */
Coord Ellipsoide::trouveIntersection(Rayon const & ray) const
{
    Coord origine( VersLocal( ray.origine() ) );
    Vecteur dir( VersLocal( ray.direction() ) );

    dir     *= VersSphere;
    origine *= VersSphere;

    double b( dir * origine );
    b *= 2.;
    double c( origine.r2() - 1. );

    MoteurResolution -> setParams( c, b, dir.norme() );

    vector< double > solutions;

    MoteurResolution -> resoud( solutions );

    vector< double >::const_iterator It, Begin( solutions.begin() ),
            End( solutions.end() );

    Coord local;

    for ( It = Begin; It != End; It++ )
    {
        if ( *It > 1e-7 )
        {
            local = origine;
            local += (*It) * dir;
            local *= VersSphere.inverse();
            return VersGlobal( local );
        }
    }

    return NullePart;
}

/** Calcule la valeur du champ de normale au point donné.
 *
 * Les composantes du vecteur normal sont simplement données par la direction
 * du vecteur \f$ \vec{OP} \f$.
 *
 * @warning Ne vérifie pas si le point est bien à la surface de l'objet.
 *
 * @param[in] P point de la surface pour lequel on désire connaître la
 * normale.
 *
 * @return la normale à l'ellipsoïde par le point \c P.
 */
VecteurU Ellipsoide::normale(Coord const & P) const
{
    return VecteurU( Centre, P );
}

/** Permet de vérifier si un point est à @c l'intérieur de l'instance.
 *
 * On passe premièrement dans le système de coordonées locales, on se réduit
 * ensuite à la sphère unité et on compare la norme du vecteur obtenu avec 1.
 *
 * @param[in] P point à tester.
 *
 * @retval true si le point est bien dans le volume, @retval false sinon.
 */
bool Ellipsoide::estDans(Coord const & P) const
{
    Vecteur local( P );

    local = VersLocal( local );

    local *= VersSphere;

    return ( local.norme() - 1. < 1.e-14  );
}

/** Permet de vérifier si un point est à la @e surface de l'instance.
 *
 * On passe premièrement dans le système de coordonées locales, on se réduit
 * ensuite à la sphère unité et on compare la norme du vecteur obtenu avec 1.
 *
 * @param[in] P point à tester.
 *
 * @retval true si le point est bien sur la surface, @retval false sinon.
 */
bool Ellipsoide::estSur(Coord const & P) const
{
    Vecteur local( P );

    local = VersLocal( local );

    local *= VersSphere;

    return ( abs( local.norme() - 1. ) < 1.e-14  );
}

/** Méthode permettant d'appliquer une texture sur l'ellipsoïde.
 *
 * On récupère les coordonnées intrinsèques, que l'on traite pour déterminer le
 * pixel à mapper.
 *
 * Si le point en question ne peut être converti en coordonnées intrinsèques,
 * alors le pixel (0,0) est utilisé.
 *
 * @param[in] pt point de la surface à mapper.
 * @param[out] i coordonnée x du pixel de l'image à plaquer au point pt.
 * @param[out] j coordonnée y du pixel de l'image à plaquer au point pt.
 * @param[in] frequence fréquence de répétition de la texture.
 * @param[in] Hauteur hauteur (en pixels) de l'image à mapper.
 * @param[in] Largeur largeur (en pixels) de l'image à mapper.
 */
void Ellipsoide::mapping(const Coord &pt, int &i, int &j,
                         const int &frequence, const int &Hauteur,
                         const int &Largeur) const
{
    Coord carte( intrinseque( pt ) );

    if ( carte == NullePart )
    {
        i = j = 0;
        return;
    }

    double phi( carte.x() / RayonX ), theta( carte.y() / RayonY ), t;

    t = double( frequence * Hauteur ) * ( theta / Forme::ValeurPi );
    j = static_cast< int >( floor( t ) );

    t = double( frequence * Largeur ) * phi / ( 2 * Forme::ValeurPi );
    i = static_cast< int >( floor( t ) );

    i = i % Largeur;
    j = j % Hauteur;
}

/** Méthode affichant les paramètres de l'instance.
 *
 * Cette méthode affiche le centre, les axes, les rayons et le volume occupé.
 *
 * @param[in] os référence sur le flot de sortie à utiliser.
 *
 * @return la référence sur le flot de sortie utilisé.
 */
ostream& Ellipsoide::affiche(ostream& os) const
{
    os << "Forme : " << nom() << endl << "Paramètres :" << endl
       << "\tCentre : " << Centre << endl
       << "\tAxeX  : " << AxeX << endl
       << "\tAxeY  : " << AxeY << endl
       << "\tAxeZ  : " << AxeZ << endl
       << "\tRayonX  : " << RayonX << endl
       << "\tRayonY  : " << RayonY << endl
       << "\tRayonZ  : " << RayonZ << endl;
    os << "Volume centré en " << PseudoCentre << ", rayon = "
       << PseudoRayon << endl;

    return os;
}

/** Méthode calculant les coordonnées intrinsèques d'un point de la surface.
 *
 * On renvoie \f$ ( r \, \phi, r \, \theta, 0 ) \f$ où \f$\phi\f$ est
 * dans l'intervalle \f$\left[ 0, 2\pi \right[\f$ et \f$\theta\f$ dans
 * \f$\left[ 0, \pi \right[\f$.
 *
 * @param[in] pt point de la surface.
 *
 * @retval NullePart si le point n'est pas sur la surface, @retval local dont
 * les coordonnées sont @f$ ( r \, \phi, r \, \theta, 0 ) @f$.
 */
Coord Ellipsoide::intrinseque(const Coord& pt) const
{
    if ( ! estSur( pt ) )
        return NullePart;

    Coord local( VersLocal( pt ) );

    local *= VersSphere;

    double phi, theta;

    theta = acos( local.z() );

    if ( abs( local.y() ) < 1.e-7 && local.x() > 1.e-7 )
        phi = 0.;
    else if ( abs( local.y() ) < 1.e-7 && local.x() < -1.e-7 )
        phi = Forme::ValeurPi;
    else
        phi = acos( local.x() / sin( theta ) );

    if ( local.y() < -1.e-7 )
        phi = 2. * Forme::ValeurPi - phi;

    return Coord( phi, theta, 0. ) * VersSphere.inverse();
}

/** Méthode calculant les coordonnées intrinsèques réduites d'un point de la
 * surface.
 *
 * On divise les deux angles par @f$ \pi @f$ pour les ramener dans l'intervalle
 * \f$\left[ 0, 2 \right[\f$ et \f$\left[ 0, 1 \right[\f$ respectivement, puis
 * on multiplie par la largeur voulue.
 *
 * @param[in] pt le point de la surface.
 * @param[in] larg
 *
 * @retval NullePart si le point n'est pas à la surface, @retval local dont les
 * coordonnées sont @f$ \frac{larg}{\pi}( r \, \phi, r \, \theta, 0 ) @f$.
 */
Coord Ellipsoide::intrinsequeReduite(const Coord& pt, const double& larg) const
{
    if ( ! estSur ( pt ) )
        return NullePart;

    Coord carte( intrinseque( pt ) );

    // Attention, theta varie seulement dans l'intervalle [0, pi[

    return carte *= ( larg / Forme::ValeurPi );
}
